/**
 * 
 */
package cn.cityhouse.avm.regression;

import java.util.ArrayList;
import java.util.Collections;

import org.apache.commons.math3.distribution.FDistribution;
import org.apache.commons.math3.distribution.NormalDistribution;
import org.apache.commons.math3.distribution.TDistribution;
import org.apache.commons.math3.exception.DimensionMismatchException;
import org.apache.commons.math3.exception.NullArgumentException;
import org.apache.commons.math3.linear.RealMatrix;
import org.apache.commons.math3.linear.RealVector;
import org.apache.commons.math3.linear.SingularValueDecomposition;
import org.apache.commons.math3.stat.regression.OLSMultipleLinearRegression;

import cn.cityhouse.avm.stat.SummaryStats;
import cn.cityhouse.avm.utils.DataFrame;
import cn.cityhouse.avm.utils.Out;

/**
 * @author <a href="mailto:wangcl@lreis.ac.cn">ChenLiang Wang</a>
 *
 */
public class AVMOLRegression extends OLSMultipleLinearRegression implements AVMAbstractRegression {

	/**
	 * 
	 */
	double[] residuals = null;
	double rsquare = Double.NaN;
	double adj_rsquare = Double.NaN;
	double d_sigma = Double.NaN;
	int int_k = -1;// the counter of paramters
	double[] d_y = null;
	double[][] d_x = null;
	double d_SSR = Double.NaN;// regression.calculateResidualSumOfSquares();
	double d_SST = Double.NaN;// regression.calculateTotalSumOfSquares();
	double d_SSE = Double.NaN;// d_SST - d_SSR;
	double d_AIC = Double.NaN;
	double d_AICc = Double.NaN;
	private String[] x_name = null;
	private String y_name = null;
	double[] d_fitted_y = null;
	double vary = Double.NaN;
	// double[] errors =
	// regression.estimateRegressionParametersStandardErrors();Returns the
	// standard errors of the regression parameters:
	double[] errors = null;
	// double[][] varb = regression.estimateRegressionParametersVariance();
	double[][] varb = null;
	// double[] betaHat = regression.estimateRegressionParameters();
	double[] betaHat = null;
	// double[] d_tstatisc_coef= t statistic (t statistic=betaHat/errors)
	double[] d_tstatisc_coef = null;
	double[] d_pv_coef = null;// Significance of Regression Coefficient ,
	// two-tailed test p value
	double d_df_coef_t = Double.NaN;
	double d_CI = 0.95;// confidence interval (CI)
	double d_fstat = 0;
	double d_pv_reg = -1;// p value of model(F test)
	double loglik = Double.NaN;
	double d_sigma_ML = Double.NaN; // d_sigma * Math.sqrt((double) (nrow -
									// int_k) / nrow);
	double[] upperBeta = null;
	double[] lowerBeta = null;

	public AVMOLRegression() {
		// TODO Auto-generated constructor stub
	}

	/**
	 * @param threshold
	 *            the compute accuracy threshold for the regression
	 */
	public AVMOLRegression(double threshold) {
		super(threshold);
		// TODO Auto-generated constructor stub
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#calculateFittedValue()
	 */
	@Override
	public RealVector calculateFittedValue() {
		// TODO Auto-generated method stub
		RealVector b = calculateBeta();
		return this.getX().operate(b);
	}
	/**
	 * <p>The condition number of a linear regression problem is a worst-case measure of
	 *  how sensitive the result is to small perturbations of the input.
	 * A large condition number (say, more than 1000) indicates the presence of
	 *  significant multicollinearity.</p>
	 * @return The condition number
	 */
	public double estimateCondNum(){
		SingularValueDecomposition svd =
	            new SingularValueDecomposition(this.getX());
		return(svd.getConditionNumber());
	}

	public double calculateLogLik() {

		int nrow = this.getNrow();
		int int_k = this.getK();
		this.d_sigma_ML = d_sigma * Math.sqrt((double) (nrow - int_k) / nrow);
		NormalDistribution p_norm = null;
		d_fitted_y = estimateFittedValue();
		if (d_fitted_y.length != nrow) {
			throw new DimensionMismatchException(d_fitted_y.length, nrow);
		}
		double loglik = 0;
		for (int i = 0; i < nrow; i++) {
			p_norm = new NormalDistribution(d_fitted_y[i], d_sigma_ML);
			loglik += p_norm.logDensity(d_y[i]);
			// System.out.println(loglik);
		}
		return loglik;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#estimateAIC()
	 */
	@Override
	public double estimateAIC() {
		// TODO Auto-generated method stub
		int nrow = this.getNrow();
		this.d_SSR = calculateResidualSumOfSquares();
		this.d_SST = calculateTotalSumOfSquares();
		this.d_SSE = d_SST - d_SSR;
		double AIC = nrow * Math.log(d_SSR / nrow) + nrow * Math.log(2.0d * Math.PI) + nrow
				+ 2.0d * (this.getK() + 1.0d);
		return AIC;
	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#estimateAICc()
	 */
	@Override
	public double estimateAICc() {
		// TODO Auto-generated method stub
		int nrow = this.getNrow();
		this.d_SSR = calculateResidualSumOfSquares();
		this.d_SST = calculateTotalSumOfSquares();
		this.d_SSE = d_SST - d_SSR;
		double d_AICc = nrow * Math.log(d_SSR / nrow) + nrow * Math.log(2.0d * Math.PI) + nrow
				+ 2.0d * nrow * (this.getK() + 1.0d) / (nrow - this.getK() - 2.0d);
		return d_AICc;
	}

	public void estimateBetaPV() {
		int nrow = this.getNrow();
		int int_k = this.getK();
		if (nrow <= 0 || int_k <= 0) {
			throw new NullPointerException();
		}
		this.d_tstatisc_coef = new double[int_k];
		this.d_pv_coef = new double[int_k];
		this.d_df_coef_t = nrow - int_k;
		TDistribution t = new TDistribution(this.d_df_coef_t);
		// double d_CI = 0.95; // confidence interval (CI)
		System.out.printf("%15s", Out.round(((1.0 - d_CI) * 50), 2) + " %:\t");
		this.lowerBeta = new double[this.betaHat.length];
		this.upperBeta = new double[this.betaHat.length];
		for (int i = 0; i < int_k; i++) {
			this.d_tstatisc_coef[i] = this.betaHat[i] / this.errors[i];
			// two-tailed test
			this.d_pv_coef[i] = 2 * (1 - t.cumulativeProbability(Math.abs(this.d_tstatisc_coef[i])));
			this.lowerBeta[i] = betaHat[i] - SummaryStats.getBetaConfidenceInterval(this.errors[i], 1 - d_CI, nrow);
			System.out.printf("|%15.6f", this.lowerBeta[i]);
		}
		System.out.println();
		System.out.printf("%15s", Out.round((100 - (1.0 - d_CI) * 50), 2) + " %:\t");
		for (int i = 0; i < int_k; i++) {
			this.upperBeta[i] = betaHat[i] + SummaryStats.getBetaConfidenceInterval(this.errors[i], 1 - d_CI, nrow);
			System.out.printf("|%15.6f", this.upperBeta[i]);
		}
		System.out.println();
		System.out.printf("%15s", "t value:\t");
		Out.printArr(this.d_tstatisc_coef);
		System.out.printf("\n%15s", "Pr(>|t|):\t");
		Out.printArr(this.d_pv_coef, "|%15.5e");
		System.out.println();

	}

	/*
	 * (non-Javadoc)
	 * 
	 * @see cn.cityhouse.avm.regression.AVMRegression2#estimateFittedValue()
	 */
	@Override
	public double[] estimateFittedValue() {
		// TODO Auto-generated method stub
		return calculateFittedValue().toArray();
	}

	public void ftestModel() {
		// double d_fstat = 0;
		int nrow = this.getNrow();
		int int_k = this.getK();
		this.d_SSR = calculateResidualSumOfSquares();
		this.d_SST = calculateTotalSumOfSquares();
		this.d_SSE = d_SST - d_SSR;
		this.d_fstat = d_SSE / d_SSR * (nrow - int_k) / (int_k - 1);
		System.out.printf("\nF-statistic: %.6f on %d and %d DF", d_fstat, (int_k - 1), (nrow - int_k));
		FDistribution f = new FDistribution(int_k - 1, nrow - int_k);
		this.d_pv_reg = 1 - f.cumulativeProbability(d_fstat);
		System.out.println(",  p-value: " + d_pv_reg);
		// System.out.println("-----\nSignif. codes: 0 '***' 0.001 '**' // 0.01
		// '*' 0.05 '.' 0.1 ' ' 1 ");
	}

	public double getD_SSR() {
		return d_SSR;
	}

	@Override
	public int getK() {
		// TODO Auto-generated method stub
		return getX().getColumnDimension();
	}

	@Override
	public int getNrow() {
		// TODO Auto-generated method stub
		return getX().getRowDimension();
	}

	@Override
	public RealMatrix getX() {
		return super.getX();
	}

	/**
	 * Return names of the independent variables
	 * 
	 * @return names of the independent variables
	 */
	public String[] getX_name() {
		ArrayList<String> variables = new ArrayList<String>();
		Collections.addAll(variables, this.x_name);
		if (!isNoIntercept()) {
			variables.add(0, "Intercept");
		}
		String[] xname = new String[x_name.length + (isNoIntercept() ? 0 : 1)];
		xname = variables.toArray(xname);
		return xname;
	}

	@Override
	public RealVector getY() {
		return super.getY();
	}

	/**
	 * Return the name of the dependent variable
	 * 
	 * @return name of the dependent variable
	 */
	public String getY_name() {
		return y_name;
	}

	/**
	 * 
	 * @param data_frame
	 *            a data object containing the variables in the model
	 * @param x_name
	 *            names of the independent variables
	 * @param y_name
	 *            name of the dependent variable
	 */
	public void loadData(DataFrame data_frame, String[] x_name, String y_name) {
		// TODO Auto-generated method stub
		if (data_frame == null || y_name == null || x_name == null) {
			throw new NullArgumentException();
		}
		double[] d_y = data_frame.getColumn(y_name);
		double[][] d_x = data_frame.getColumns(x_name);// subset
		newSampleData(d_y, d_x);
		this.x_name = x_name;
		this.y_name = y_name;
	}

	/**
	 * @param d_y
	 *            double array of the dependent variable
	 * @param d_x
	 *            double array of the independent variables
	 * @param x_name
	 *            names of the independent variables
	 * @param y_name
	 *            name of the dependent variable
	 */
	public void loadData(double[] d_y, double[][] d_x, String[] x_name, String y_name) {
		// TODO Auto-generated method stub
		if (d_y == null || d_x == null || x_name == null || y_name == null) {
			throw new NullArgumentException();
		}
		newSampleData(d_y, d_x);
		this.x_name = x_name;
		this.y_name = y_name;
	}

	@Override
	public void printSummary() {
		// TODO Auto-generated method stub
		betaHat = estimateRegressionParameters();
		residuals = estimateResiduals();
		System.out.println("Estimates the residuals, ie u=y-X*b: ");
		// Out.printArr(residuals, df);
		SummaryStats.summaryStats(residuals);
		vary = estimateRegressandVariance();
		// System.out.print("Returns the variance of the regression Var(y): ");
		// System.out.println(vary);
		errors = estimateRegressionParametersStandardErrors();
		varb = estimateRegressionParametersVariance();
		rsquare = calculateRSquared();
		adj_rsquare = calculateAdjustedRSquared();
		d_sigma = estimateRegressionStandardError();
		System.out.printf("\n%15s", "Coefficients:\t");
		Out.printArr(getX_name());
		System.out.printf("\n%15s", "Estimate:\t");
		Out.printArr(betaHat);
		System.out.printf("\n%15s", "Std. Error:\t");
		Out.printArr(errors);
		this.d_fitted_y = this.estimateFittedValue();
		System.out.println();
		this.estimateBetaPV();
		this.ftestModel();
		System.out.printf("Residual sum of squares:\t%15.6f\n", this.d_SSR);
		System.out.println("Residual standard error: " + d_sigma);
		System.out.printf("AIC: %.6f , AICc: %.6f , BIC: %.6f\n", 
				this.estimateAIC(), this.estimateAICc(),this.estimateBIC());
		System.out.print("Multiple R-squared: " + rsquare);
		System.out.println("\tAdjusted R-squared: " + adj_rsquare);
	}

	@Override
	public double estimateBIC() {
		// TODO Auto-generated method stub
		int nrow = this.getNrow();
		this.d_SSR = calculateResidualSumOfSquares();
		double BIC = nrow * Math.log(d_SSR / nrow) + nrow * Math.log(2.0d * Math.PI) + nrow
				+ Math.log(nrow) * (this.getK() + 1.0d);		
		//BIC=-2.0d*this.calculateLogLik()+Math.log(nrow)*(this.getK()+1.0d);
		//BIC=this.estimateAIC()-2.0d*(this.getK()+1.0d)+Math.log(nrow) * (this.getK() + 1.0d);
		return BIC;
	}

}
